RRST : Analysis - sample excluding outlier

Author

Natascha Stoffel

Published

November 22, 2024

0) Data and Preparation

Respiratory Interoception; Analysis of Site Switzerland

Analysis of participants of site CH from 15.10.23 until 10.09.2024;

# remove everything that we stored in the environment
rm(list=ls()) 

# Load necessary libraries
library(readr)
library(dplyr)
library(rstatix)
library(gtsummary)
library(ggplot2)
library(ggpubr)
library(openxlsx)
library(gt)
library(stats)
library(grDevices)
library(car)
library(stringr)  
library(lubridate)
library(Hmisc)

### 0.0 PREPARATION -----------------------------------------------------------
# load prepared and cleaned file from 00_preparation: df_yyyymmdd.csv
setwd("/Users/nataschastoffel/Documents/GitHub/interoception_NS")
df_RRST <- read_csv("/Users/nataschastoffel/Documents/GitHub/interoception_NS/data/processed/11-2024/Data.RRST_FINAL_20241120.csv")
df <- df_RRST[-c(1)] # remove additional column 

# excluded_pcodes <- c(14)  # pcodes that we want to exclude due to our exclusion criteria (that they were not able to perform the task adequetely; due to not understanding/performing with a 0.0000 sensitivity)
# Filter out the rows with the specified pcodes
# df <- df %>% filter(!pcode %in% excluded_pcodes)

#Set path for output files; working directory to the correct folder
setwd("~/Documents/GitHub/interoception_NS/output/RRST/11-2024")

df <- df %>%
  mutate(anx_dep_SUM = stai_t + bdi)# create one variable for the trait anxiety and depression to use as one variable in the model


# factor variables
factor_variables <- c('group', 'sex', 'psychotropic_medication', 'menopause', 'contraception')
# Apply as.numeric to the selected columns
df[factor_variables] <- lapply(df[factor_variables], as.factor)

# numeric variables
numeric_variables <- c('age', 'bmi', 'bdi', 'stai_s', 'stai_t', 'maia_total', 'anx_dep_SUM',
                     'sdq', 'ctq_total', 'ias', 'rrst_sensitivity', 
                     'traitPE_maiaRRST', 'traitPE_iasRRST', 'metascore')
# Apply as.numeric to the selected columns
df[numeric_variables] <- lapply(df[numeric_variables], as.numeric)

# define required subgroups
df_female <- df[df$sex == "1", ]
df_male <- df[df$sex == "2", ]

df_HC <- df[df$group == "0", ]
df_FND <- df[df$group == "1", ]

Identifying & Argueing for the Outlier to exclude

# Testing Outlier using boxplot (standard deviation)
df %>% 
  group_by(group) %>%
  identify_outliers(rrst_sensitivity)
# A tibble: 1 × 77
  group pcode date         age sex   handedness diagnosis medication.x smoke_crf
  <fct> <dbl> <date>     <dbl> <fct>      <dbl>     <dbl>        <dbl>     <dbl>
1 1        63 2024-05-08    38 1              1         5            1         0
# ℹ 68 more variables: drugs_crf <dbl>, bmi <dbl>, bdi <dbl>, stai_t <dbl>,
#   stai_s <dbl>, ias <dbl>, maia_note <dbl>, maia_distr <dbl>,
#   maia_worry <dbl>, maia_attreg <dbl>, maia_aware <dbl>, maia_sfreg <dbl>,
#   maia_body <dbl>, maia_trust <dbl>, maia_total <dbl>, sdq <dbl>,
#   contraception <fct>, menopause <fct>, mc_phase <dbl>, sfmdrs <dbl>,
#   seiz <dbl>, cgi <dbl>, ctq_emoab <dbl>, ctq_physab <dbl>, ctq_sexab <dbl>,
#   ctq_emoneg <dbl>, ctq_physneg <dbl>, ctq_minim <dbl>, ctq_total <dbl>, …
# one outlier identified = P063 for rrst_sensitivity 

#### Box-Plot with Jitter
plotBoxplotGroups <- function(df, var , tit ){ 
  plot <- ggplot(df, aes(x = group, y = var, fill = group)) + 
    geom_boxplot(show.legend = TRUE) + # Show legend for boxplot
    labs(y = "Score") + 
    theme_classic() + 
    scale_fill_manual(values = c("#868686FF", "#BB4038")) + 
    scale_x_discrete(labels = c("HC", "FND")) +  
    guides(fill = "none") + # Remove the legend
    ggtitle(tit) +
    xlab("") +
    geom_jitter(shape = 16, position = position_jitter(0.2))
  
  return(plot)
}

p1 <- plotBoxplotGroups(df, (df$rrst_sensitivity), "Respiration Sensitivity") + 
  stat_compare_means(method = "t.test", paired = FALSE, label.x = 1.3)
p1

# Filter the rows where rrst_sensitivity is greater than 8 and display pcode (visually detected as extreme)
filtered_pcodes <- df[df$rrst_sensitivity > 8, "pcode"]
print(filtered_pcodes) 
# A tibble: 1 × 1
  pcode
  <dbl>
1    63
# P063

# Calculate mean and standard deviation of rrst_sensitivity
mean_rrst <- mean(df$rrst_sensitivity)
sd_rrst <- sd(df$rrst_sensitivity)
# Define the threshold for filtering (mean + 2.5 * SD)
threshold <- mean_rrst + 2.5 * sd_rrst # 7.07
# Filter the pcodes for rrst_sensitivity greater than the threshold
filtered_pcodes <- df[df$rrst_sensitivity > threshold, "pcode"]
# View the result
filtered_pcodes
# A tibble: 1 × 1
  pcode
  <dbl>
1    63
# P063

# Testing the Assumptions in the main model (rrst sensitivity predicted by clinical variables) 
# model:
lmrrst <- lm(
  rrst_sensitivity ~ group + sex + bmi + psychotropic_medication + anx_dep_SUM + ctq_total + sdq,
  data = df)
# testing assumptions
plot(lmrrst, which=1) # testing the assumptions of Linearity and Homoscedasticity with residuals vs fitted plot: the line should more or less be straight: observation 56 is a clear outlier

plot(lmrrst, which=2) # testing the assumptions normality of residuals with  standardized quantiles: the data points should more or less on the line: observation 56 is a clear outlier (not as extreme but also noted: 81, 65)

plot(lmrrst, which=3) # testing the assumption of homoscedasticity with scale-location plot--> if there is an incline then we should consider log transformation: observation 56 is a clear outlier (not as extreme but also noted: 81, 65)

plot(lmrrst, which=5) # testing leverage and influence = when an outlier has a lot of leverage over the regression line (high leverage & large residual; these points will lie close to the right top or bottom corners): observation 56 is a clear outlier (not as extreme but also noted: 81, 6)

# observation 56 (consistent outlier) = P063

df_exclusion <- df[df$pcode != "63", ]

p1b <- plotBoxplotGroups(df_exclusion, (df_exclusion$rrst_sensitivity), "Respiration Sensitivity") + 
  stat_compare_means(method = "t.test", paired = FALSE, label.x = 1.3)
print(p1b)

p1b

# Summary of Change - depending on Exclusion of Outlier
summary(df$rrst_sensitivity)  # no sig difference between group (cut of 0.05)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   2.458   3.280   3.282   4.350   8.752 
summary(df_exclusion$rrst_sensitivity) # sig difference between group (cut of 0.05)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   2.429   3.277   3.222   4.307   6.242 
# for statistical reasons we EXCLUDE P063 from the dataset
df <- df[df$pcode != "63", ]

1) Having a look at the Data

# Normality of Data -------------------
shapiro.test(df$rrst_sensitivity) #  normally distributed
shapiro.test(df$metascore)# NOT normally distributed
shapiro.test(df$metascore[df$group == "1"]) # patients within themselves are NOT normally distributed
shapiro.test(df$metascore[df$group == "0"]) # controls within themselves are normally distributed

shapiro.test(df$age) # normally distributed

shapiro.test(df$bdi) #  NOT normally distributed : p < 0.05
shapiro.test(df$bdi[df$group == "1"]) # patients within themselves are normally distributed
shapiro.test(df$bdi[df$group == "0"]) # controls within themselves are NOT normally distributed

shapiro.test(df$stai_t) #  normally distributed 

shapiro.test(df$stai_s) #  NOT normally distributed 
shapiro.test(df$stai_s[df$group == "1"]) # patients within themselves are normally distributed
shapiro.test(df$stai_s[df$group == "0"]) # controls within themselves are NOT normally distributed : p < 0.05

shapiro.test(df$ctq_total) #  NOT normally distributed : p < 0.05
shapiro.test(df$ctq_total[df$group == "1"]) # patients within themselves are also NOT normally distributed
shapiro.test(df$ctq_total[df$group == "0"]) # controls within themselves are also NOT normally distributed

shapiro.test(df$sdq) #  NOT normally distributed : p < 0.05
shapiro.test(df$sdq[df$group == "1"]) # patients within themselves are NOT normally distributed
shapiro.test(df$sdq[df$group == "0"]) # controls within themselves are NOT normally distributed : p < 0.05

shapiro.test(df$maia_total) # normally distributed 
shapiro.test(df$traitPE_maiaRRST) # normally distributed with p < 0.05

shapiro.test(df$ias) # NOT normally distributed 
shapiro.test(df$traitPE_iasRRST) # normally distributed with p < 0.05

shapiro.test(df$bmi) # NOT normally distributed with p < 0.05
shapiro.test(df$bmi[df$group == "1"]) # within patients bmi is normally distributed 
shapiro.test(df$bmi[df$group == "0"]) # within control bmi is NOT normally distributed with p < 0.05


## 1.1 testing co-dependency of variables by creating a correlation matrix -----
# in order to run a correlation, all variables must be numeric (but we create new variables, so it wont mess with our variables that we defiend as factors)
corr_matrix <- cor(df[, c("age","bmi", "bdi", "stai_s", "stai_t", "ctq_total",
                          "ias", "maia_total", "sdq", "rrst_sensitivity", 
                          "traitPE_iasRRST", "traitPE_maiaRRST", "metascore")], use = "complete.obs")

# Set the cutoff for high correlations, that we then want to exclude
cutoff <- 0.7

# Find the indices of correlations that are greater than 0.8 and exclude self-correlations (diagonal)
high_corr <- which(abs(corr_matrix) > cutoff & abs(corr_matrix) < 1, arr.ind = TRUE)

# Display the pairs of variables with high correlations
high_corr_pairs <- data.frame(
  Var1 = rownames(corr_matrix)[high_corr[, 1]],
  Var2 = colnames(corr_matrix)[high_corr[, 2]],
  Correlation = corr_matrix[high_corr]
)

# View the pairs
print(high_corr_pairs) # stai and bdi are highly correlated, which we do accept though as this is a common clinical happening
cor.test(df$bdi, df$stai_t) # reason for using sum scores

### 1.2 Group Differences as Baseline of Variables of Interst ---------------------------------------------
# parametric t-test for NORMALLY distributed data
res<-t.test(df$rrst_sensitivity[df$group=="1"], df$rrst_sensitivity[df$group=="0"])
res # SIG difference between groups
res<-t.test(df$stai_t[df$group=="1"], df$stai_t[df$group=="0"])
res # SIG difference between groups with p < 0.05
res<-t.test(df$maia_total[df$group=="1"], df$maia_total[df$group=="0"])
res # SIG difference between groups with p < 0.05
res<-t.test(df$traitPE_iasRRST[df$group=="1"], df$traitPE_iasRRST[df$group=="0"])
res # no difference between groups with p < 0.05
res<-t.test(df$traitPE_maiaRRST[df$group=="1"], df$traitPE_maiaRRST[df$group=="0"])
res # no sig difference between groups

# non paramatric wilcox test for NOT normally distributed data 
res<-wilcox.test(df$metascore[df$group=="1"], df$metascore[df$group=="0"])
res # no sig difference between groups
res<-wilcox.test(df$ias[df$group=="1"], df$ias[df$group=="0"])
res # SIG difference between groups
res<-wilcox.test(df$bdi[df$group=="1"], df$bdi[df$group=="0"])
res # SIG difference between groups with p < 0.05
res<-wilcox.test(df$stai_s[df$group=="1"], df$stai_s[df$group=="0"])
res # SIG difference between groups with p < 0.05
res<-wilcox.test(df$ctq_total[df$group=="1"], df$ctq_total[df$group=="0"])
res # SIG difference between groups with p < 0.05
res<-wilcox.test(df$sdq[df$group=="1"], df$sdq[df$group=="0"])
res # SIG difference between groups with p < 0.05
res<-wilcox.test(df$ias[df$group=="1"], df$ias[df$group=="0"])
res # no sig difference between groups
res<-wilcox.test(df$bmi[df$group=="1"], df$bmi[df$group=="0"])
res # no sig difference between groups

Demographic Summary Table

# overview across group and sex distribution
tbl <- df %>%
  mutate(sex = factor(sex, levels = c("1", "2"), labels = c("female", "male"))) %>%
  mutate(group = factor(group, levels = c("0", "1"), labels = c("HC", "FND")))
tbl1<-table(tbl$group,tbl$sex)
tbl1
     
      female male
  HC      35   13
  FND     31   12
# overview for mean age across group
summary_stats <- tbl %>%
  group_by(group) %>%
  get_summary_stats(age, type = "full") # 'full' includes mean, SD, min, max, etc.
print(summary_stats)
# A tibble: 2 × 14
  group variable     n   min   max median    q1    q3   iqr   mad  mean    sd
  <fct> <fct>    <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 HC    age         48    18    67     37  28      50  22    14.1  37.8  13.0
2 FND   age         43    18    64     39  31.5    47  15.5  11.9  38.5  11.8
# ℹ 2 more variables: se <dbl>, ci <dbl>
res<-t.test(df$age[df$group == "0"], df$age[df$group == "1"]) # 
res # no sig difference between the groups in age

    Welch Two Sample t-test

data:  df$age[df$group == "0"] and df$age[df$group == "1"]
t = -0.25959, df = 88.975, p-value = 0.7958
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -5.849281  4.497537
sample estimates:
mean of x mean of y 
 37.81250  38.48837 
# SUMMARY TABLE DEMOGRAPHICS -----
df_tbl <- df %>%
  mutate(sex = factor(sex, levels = c("1", "2"), labels = c("male", "female"))) %>%
  mutate(group = factor(group, levels = c("0", "1"), labels = c("HC", "FND"))) %>%
  mutate(psychotropic_medication = factor(psychotropic_medication, levels = c("0", "1"), labels = c("no", "yes"))) 

library(gtsummary)
library(dplyr)

# Select relevant columns from df_tbl
Data.t <- select(df_tbl, c("group", "sex", "age", "psychotropic_medication", "bmi", "bdi", "stai_t", "ctq_total", "sdq"))

# Create the summary table
tbl1 <- Data.t %>%
    tbl_summary(
    by = group,
    missing = "no",
    statistic = list(
      age = "{mean} ({sd})",
      sex = "{n} ({p})",
      psychotropic_medication = "{n} ({p})",
      bmi = "{median} ({p25}, {p75})",
      bdi = "{median} ({p25}, {p75})",
      stai_t = "{mean} ({sd})",
      ctq_total = "{median} ({p25}, {p75})",
      sdq = "{median} ({p25}, {p75})"
    ),
    digits = list(age = 1, sex = 1, psychotropic_medication = 1, bmi = 1, bdi = 1, stai_t=1, ctq_total=1, sdq=1),
    label = list(age = "Age, mean (SD)",
                 sex = "Sex, count (%)",
                 psychotropic_medication= "Intake of psychotropic Medication, count (%)",
                 bmi = "Body Mass Index kg/m², median (IQR)",
                 bdi = "Depression using BDI-II, median (IQR), ",
                 stai_t = "Trait Anxiety using STAI, median (IQR)",
                 ctq_total = "Childhood Trauma, using CTQ total, median (IQR)",
                 sdq = "Dissociation using SDQ-20, median (IQR)")
  ) %>%
  modify_header(list(label ~ "**Variable**")) %>%
  modify_caption("**Demographics Overview**") %>%
  modify_spanning_header(c("stat_1", "stat_2") ~ "**Group**") %>%
  add_p( test = list(
        age ~ "t.test",            # t-test for numeric variables (even though NOT normally distributed; rather use wilcox and then the median instead of the mean??)
        bmi ~ "wilcox.test",            # t-test for bmi
        sex ~ "chisq.test",        # chi-square test for factor variables
        psychotropic_medication ~ "chisq.test",  # Fisher's exact test for small samples
        bdi ~ "wilcox.test",
        stai_t ~ "t.test",
        ctq_total = "wilcox.test",
        sdq="wilcox.test"
      )) %>% #add p-value# add p-value (in case different variable should be tested differently, create a list and specify like that: add_p(test = list(age = "t.test", bmi = "wilcox.test"))
  add_overall(last = FALSE)

# Display the table
tbl1
Demographics Overview
Variable Overall, N = 911 Group p-value2
HC, N = 481 FND, N = 431
Sex, count (%)


>0.9
    male 66.0 (72.5) 35.0 (72.9) 31.0 (72.1)
    female 25.0 (27.5) 13.0 (27.1) 12.0 (27.9)
Age, mean (SD) 38.1 (12.4) 37.8 (13.0) 38.5 (11.8) 0.8
Intake of psychotropic Medication, count (%) 21.0 (23.6) 2.0 (4.3) 19.0 (44.2) <0.001
Body Mass Index kg/m², median (IQR) 23.7 (21.7, 26.3) 22.7 (21.7, 25.1) 25.3 (21.8, 27.8) 0.019
Depression using BDI-II, median (IQR), 8.0 (3.0, 15.0) 4.0 (1.0, 8.0) 15.0 (9.5, 23.5) <0.001
Trait Anxiety using STAI, median (IQR) 41.6 (12.1) 37.3 (10.3) 46.5 (12.1) <0.001
Childhood Trauma, using CTQ total, median (IQR) 37.0 (30.0, 54.5) 34.0 (29.0, 41.3) 46.0 (35.0, 60.5) 0.001
Dissociation using SDQ-20, median (IQR) 27.0 (22.0, 36.0) 23.0 (21.0, 26.0) 36.0 (28.5, 45.0) <0.001
1 n (%); Mean (SD); Median (IQR)
2 Pearson’s Chi-squared test; Welch Two Sample t-test; Wilcoxon rank sum test

FND Summary Table

# 1.4 SUMMARY TABLE for FND symptoms -----
# create variable of duation of months
df_FND <- df_FND %>%
  mutate(duration_months = interval(date_symptom_onset, date) %/% months(1)) # calculate month duration

as.numeric(df_FND$duration_months)
 [1]  32 168  39  26  43  41  25  28  25 305  72  15  12  37  36  71  69  63 173
[20]  59 207 139  71  52 259  34  39 113  34  60  22  20  41  58  31  29  39 161
[39]  89  74 124  11  69  29
# Select only the defined variables
df_tbl <- df_FND %>%
  mutate(sex = factor(sex, levels = c("1", "2"), labels = c("female", "male"))) 

library(gtsummary)
library(dplyr)

# Select relevant columns from df_tbl
Data.t <- select(df_tbl, c("sex", "duration_months", "fds", "motor", "weakness", "sensory", "pppd", "cognitive"))

# Create the summary table
tbl <- Data.t %>%
  tbl_summary(
    by = sex,
    missing = "no",
    statistic = list(
      duration_months = "{mean} ({sd})",
      fds = "{n} ({p})",
      motor = "{n} ({p})",
      weakness = "{n} ({p})",
      sensory = "{n} ({p})",
      pppd = "{n} ({p})",
      cognitive = "{n} ({p})"
      ),
    digits = list(duration_months = 1,fds = 1,motor = 1,weakness = 1,sensory = 1,pppd = 1,cognitive = 1),
    label = list(
      sex = "Sex",
      duration_months = "Symptom Duration [in months]",
      fds = "Functional Dissociative Seizures [yes]",
      motor = "motor + symptoms [yes]",
      weakness = "motor - symptoms [yes]",
      sensory = "sensory symptoms [yes]",
      pppd = "dizziness (PPPD) [yes]",
      cognitive = "cognitive symptoms [yes]"
    )
  ) %>%
  modify_header(list(label ~ "**Variable**")) %>%
  modify_caption("**Clinical Characteristics**") %>%
  modify_spanning_header(c("stat_1", "stat_2") ~ "**Group**") %>%
  add_overall(last = FALSE)

# Print the summary table
tbl
Clinical Characteristics
Variable Overall, N = 441 Group
female, N = 321 male, N = 121
Symptom Duration [in months] 71.5 (66.1) 75.6 (71.1) 60.3 (51.3)
Functional Dissociative Seizures [yes] 8.0 (18.2) 8.0 (25.0) 0.0 (0.0)
motor + symptoms [yes] 20.0 (45.5) 13.0 (40.6) 7.0 (58.3)
motor - symptoms [yes] 30.0 (68.2) 21.0 (65.6) 9.0 (75.0)
sensory symptoms [yes] 20.0 (45.5) 16.0 (50.0) 4.0 (33.3)
dizziness (PPPD) [yes] 2.0 (4.5) 2.0 (6.3) 0.0 (0.0)
cognitive symptoms [yes] 4.0 (9.1) 4.0 (12.5) 0.0 (0.0)
1 Mean (SD); n (%)

2) Group Differences - Baseline

Here we visualize the different distribution per group across all variables

df$group <- as.factor(df$group)

plotViolinGroups <- function(df, var , tit ){ 
  plot <-ggplot(df, aes(x=group, y=var, fill=group)) + # this is just the ground structure. We have group on x axis, our variable score on y, and we want to seperate the groups. 
    geom_violin(show.legend = FALSE, width = 0.7, trim = FALSE, alpha = 0.8)+ # here we say what kind of plot we want to make; namely violin to show the distrubition of the original data (width to make it less wide and thus prettier) 
    stat_summary(show.legend = FALSE, fun = median, geom = "point", color = "black", size = 2) +  # Add median point without legend
    labs(y = "Score") + # These are just the different titles
    theme_classic()+ # here we have the background, simple white :)
    scale_fill_manual(values = c("#868686FF", "#BB4038"))+ # google HEX codes and find your own colours. 
    scale_x_discrete(labels=c("Controls", "Patients"))+  # these are our labels for the x axis
    guides(fill=guide_legend(title="Timepoint"))+ # and a title for our legend
    ggtitle(tit) +
    xlab("")
  
  return(plot)
}

p1 <- plotViolinGroups(df,df$bdi,"BDI - Depression")+ ggpubr::stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
p1

p2 <- plotViolinGroups(df,df$stai_s,"STAI - State Anxiety")+ stat_compare_means(method = "t.test", paired = FALSE, label.x = 1.3)
p3 <- plotViolinGroups(df,df$stai_t,"STAI - Trait Anxiety")+ stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
gridExtra::grid.arrange(p2, p3, ncol = 2)

p4 <- plotViolinGroups(df,df$ctq_total,"CTQ - Childhood Trauma")+ stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
p4

p5 <- plotViolinGroups(df,df$sdq,"SDQ-20 - Dissociation")+ stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
p5

p6 <- plotViolinGroups(df,df$maia_total,"MAIA - Interoceptive Awareness")+ stat_compare_means(method = "t.test", paired = FALSE, label.x = 1.3)
p6

p7 <- plotViolinGroups(df,df$ias,"IAS - Interoceptive Accuracy")+ stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
p7

p8 <- plotViolinGroups(df,df$traitPE_maiaRRST,"TraitPE MAIA-RRST")+ stat_compare_means(method = "t.test", paired = FALSE, label.x = 1.3)
p9 <- plotViolinGroups(df,df$traitPE_iasRRST,"TraitPE IAS-RRST")+ stat_compare_means(method = "t.test", paired = FALSE, label.x = 1.3)
gridExtra::grid.arrange(p8, p9, ncol = 2)

p10 <- plotViolinGroups(df,df$rrst_sensitivity,"Respiratory Sensitivity")+ stat_compare_means(method = "t.test", paired = FALSE, label.x = 1.3)
p11 <- plotViolinGroups(df,df$metascore,"Respiratory Metacognition")+ stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
gridExtra::grid.arrange(p10, p11, ncol = 2)

# Table for dependent Variables: group differences
df_tbl <- df %>%
  mutate(group = factor(group, levels = c("0", "1"), labels = c("HC", "FND"))) 
  
Data.t <- select(df_tbl, c("group", "maia_total","traitPE_maiaRRST", "ias", "traitPE_iasRRST", "rrst_sensitivity" , "metascore"))
library(gtsummary)
tbl2 <- Data.t %>%
  tbl_summary(
    by = group,
    missing = "no",
    statistic = list(
      maia_total = "{mean} ({sd})",
      traitPE_maiaRRST = "{mean} ({sd})",
      ias = "{median} ({IQR})",
      traitPE_iasRRST = "{mean} ({sd})",
      rrst_sensitivity = "{mean} ({sd})",
      metascore = "{median} ({IQR})"
    ),
    digits = list(maia_total = 1, traitPE_maiaRRST = 1, ias = 1, traitPE_iasRRST = 1, 
                  rrst_sensitivity = 2, metascore = 2),
    label = list(maia_total= "Interoceptive Awareness Self-Report - MAIA, Mean, (SD)",
                 traitPE_maiaRRST= "Respiratory Trait-Prediction-Error, MAIA-RRST, Mean (SD)",
                 ias= "Interoceptive Accuracy Self-Report - IAS, Median (IQR)",
                 traitPE_iasRRST = "Respiratory Trait-Prediction-Error, IAS-RRST, Mean (SD)",
                 rrst_sensitivity= "Respiratory Sensitivity - RRST, Mean (SD)",
                 metascore= "Respiratory Metacognition - AUC, Median (IQR)"
                 )
  ) %>%
  modify_header(list(label ~ "**Variable**")) %>%
  modify_caption("**Interoception Overview**") %>%
  modify_spanning_header(c("stat_1", "stat_2") ~ "**Group**") %>%
  add_p(
    test = list(
      maia_total = "t.test",            
      traitPE_maiaRRST = "t.test",            
      ias = "wilcox.test",        
      traitPE_iasRRST = "t.test", 
      rrst_sensitivity = "t.test",
      metascore = "wilcox.test")
    ) %>%
    add_overall(last = FALSE)
tbl2
Interoception Overview
Variable Overall, N = 911 Group p-value2
HC, N = 481 FND, N = 431
Interoceptive Awareness Self-Report - MAIA, Mean, (SD) 22.1 (5.8) 24.1 (4.5) 19.9 (6.3) <0.001
Respiratory Trait-Prediction-Error, MAIA-RRST, Mean (SD) 0.0 (1.2) 0.2 (1.1) -0.1 (1.4) 0.2
Interoceptive Accuracy Self-Report - IAS, Median (IQR) 85.0 (18.5) 89.0 (14.3) 80.0 (25.5) 0.018
Respiratory Trait-Prediction-Error, IAS-RRST, Mean (SD) 0.0 (1.2) 0.1 (0.9) -0.1 (1.4) 0.4
Respiratory Sensitivity - RRST, Mean (SD) 3.22 (1.41) 3.53 (1.13) 2.88 (1.62) 0.032
Respiratory Metacognition - AUC, Median (IQR) 0.57 (0.06) 0.57 (0.05) 0.59 (0.08) 0.9
1 Mean (SD); Median (IQR)
2 Welch Two Sample t-test; Wilcoxon rank sum test; Wilcoxon rank sum exact test

3) Interoception

Here we model interoceptive variables using group and clinical predictors

# predicting group with interoceptive variables and most important covariates
glm <- glm(formula = group ~  sex + rrst_sensitivity  + metascore + maia_total + ias + anx_dep_SUM + ctq_total + sdq,
            data = df,
            family = binomial) # as group has two levels
summary(glm) # SDQ predict group

Call:
glm(formula = group ~ sex + rrst_sensitivity + metascore + maia_total + 
    ias + anx_dep_SUM + ctq_total + sdq, family = binomial, data = df)

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)   
(Intercept)      -11.96648    5.20705  -2.298  0.02155 * 
sex2               1.07654    0.80022   1.345  0.17852   
rrst_sensitivity  -0.12087    0.28041  -0.431  0.66642   
metascore         11.45437    6.90159   1.660  0.09698 . 
maia_total        -0.06882    0.07490  -0.919  0.35821   
ias               -0.01992    0.02248  -0.886  0.37556   
anx_dep_SUM        0.02600    0.02280   1.140  0.25420   
ctq_total          0.02127    0.02130   0.998  0.31806   
sdq                0.21450    0.06744   3.181  0.00147 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 125.878  on 90  degrees of freedom
Residual deviance:  64.113  on 82  degrees of freedom
AIC: 82.113

Number of Fisher Scoring iterations: 6
# testing for best model fit with covariates at hand
lmrrst1 <- lm(rrst_sensitivity ~ group + anx_dep_SUM + psychotropic_medication, data = df)
lmrrst2 <- lm(rrst_sensitivity ~ group + anx_dep_SUM + ctq_total + sdq, data = df)
lmrrst3 <- lm(rrst_sensitivity ~ group + sex + psychotropic_medication + anx_dep_SUM + ctq_total + sdq, data = df)
lmrrst4 <- lm(rrst_sensitivity ~ group + bmi + psychotropic_medication + anx_dep_SUM + ctq_total + sdq, data = df)
lmrrst5 <- lm(rrst_sensitivity ~ group + sex + bmi + psychotropic_medication + anx_dep_SUM + ctq_total, data = df)
lmrrst6 <- lm(rrst_sensitivity ~ group + sex + bmi + psychotropic_medication + anx_dep_SUM + ctq_total + sdq, data = df)
lmrrst7 <- lm(rrst_sensitivity ~ group + bmi + psychotropic_medication + anx_dep_SUM + ctq_total + sdq, data = df)
lmrrst8 <- lm(rrst_sensitivity ~ group + sex + age + bmi + psychotropic_medication + anx_dep_SUM + ctq_total + sdq, data = df)

AIC(lmrrst1, lmrrst2, lmrrst3, lmrrst4, lmrrst5, lmrrst6, lmrrst7, lmrrst8) # model 3, 6 and 8 are low
        df      AIC
lmrrst1  5 308.6969
lmrrst2  6 316.5310
lmrrst3  8 299.3683
lmrrst4  8 305.6340
lmrrst5  8 307.8920
lmrrst6  9 300.6565
lmrrst7  8 305.6340
lmrrst8 10 299.0372
anova(lmrrst3, lmrrst6, lmrrst8) # testing whether it adds sig value if we change from 3 to 6 to 8; no--> model 3 chosen
Analysis of Variance Table

Model 1: rrst_sensitivity ~ group + sex + psychotropic_medication + anx_dep_SUM + 
    ctq_total + sdq
Model 2: rrst_sensitivity ~ group + sex + bmi + psychotropic_medication + 
    anx_dep_SUM + ctq_total + sdq
Model 3: rrst_sensitivity ~ group + sex + age + bmi + psychotropic_medication + 
    anx_dep_SUM + ctq_total + sdq
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1     82 125.80                              
2     81 124.80  1    1.0022 0.6691 0.41579  
3     80 119.82  1    4.9731 3.3203 0.07217 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# predciting rrst sensitvity with covariates
lm<-lm(rrst_sensitivity ~ group + sex + psychotropic_medication + anx_dep_SUM + ctq_total + sdq,
  data = df)
summary(lm) # SEX, MEDICATION, CTQ and SDQ sig associated with RRST-sensitivity

Call:
lm(formula = rrst_sensitivity ~ group + sex + psychotropic_medication + 
    anx_dep_SUM + ctq_total + sdq, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.3161 -0.7655  0.0192  0.7099  3.1040 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)               3.183104   0.521704   6.101 3.32e-08 ***
group1                   -0.136712   0.349005  -0.392  0.69628    
sex2                      0.745765   0.303654   2.456  0.01616 *  
psychotropic_medication1 -1.068365   0.400520  -2.667  0.00921 ** 
anx_dep_SUM               0.012291   0.008244   1.491  0.13982    
ctq_total                 0.018289   0.008033   2.277  0.02541 *  
sdq                      -0.040401   0.013959  -2.894  0.00487 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.239 on 82 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.2914,    Adjusted R-squared:  0.2396 
F-statistic: 5.621 on 6 and 82 DF,  p-value: 6.316e-05
# RRST sensitivity and separate CTQ subscores
lmrrst <- lm(
  rrst_sensitivity ~ group + ctq_physab + ctq_emoab + ctq_sexab + ctq_emoneg + ctq_physneg,
  data = df)
summary(lmrrst)  # no subscale in intself associated with RRST

Call:
lm(formula = rrst_sensitivity ~ group + ctq_physab + ctq_emoab + 
    ctq_sexab + ctq_emoneg + ctq_physneg, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.84069 -0.82616 -0.02228  0.87975  3.15799 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.33514    0.42522   7.843 1.25e-11 ***
group1      -0.73021    0.31848  -2.293   0.0244 *  
ctq_physab   0.04274    0.06367   0.671   0.5039    
ctq_emoab   -0.02974    0.05494  -0.541   0.5897    
ctq_sexab    0.01818    0.04504   0.404   0.6874    
ctq_emoneg   0.01242    0.04636   0.268   0.7895    
ctq_physneg -0.00417    0.06008  -0.069   0.9448    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.414 on 84 degrees of freedom
Multiple R-squared:  0.06378,   Adjusted R-squared:  -0.003093 
F-statistic: 0.9537 on 6 and 84 DF,  p-value: 0.4616
# predciting  rrst metascore with covariates
lm<-lm(metascore ~ group + sex + psychotropic_medication + anx_dep_SUM + ctq_total + sdq,
  data = df)
summary(lm) #  SDQ and CTQ sig associated with METASCORE

Call:
lm(formula = metascore ~ group + sex + psychotropic_medication + 
    anx_dep_SUM + ctq_total + sdq, data = df)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.211016 -0.025141  0.004379  0.035825  0.103220 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)               5.910e-01  2.204e-02  26.817  < 2e-16 ***
group1                    1.805e-02  1.474e-02   1.224  0.22438    
sex2                     -1.013e-02  1.283e-02  -0.789  0.43217    
psychotropic_medication1 -2.302e-02  1.692e-02  -1.361  0.17737    
anx_dep_SUM               8.034e-05  3.482e-04   0.231  0.81812    
ctq_total                 6.564e-04  3.393e-04   1.934  0.05651 .  
sdq                      -1.620e-03  5.897e-04  -2.747  0.00739 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.05232 on 82 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.1573,    Adjusted R-squared:  0.09559 
F-statistic:  2.55 on 6 and 82 DF,  p-value: 0.02582
# RRST metascore and separate CTQ subscores
lm <- lm(
  metascore ~ group + ctq_physab + ctq_emoab + ctq_sexab + ctq_emoneg + ctq_physneg,
  data = df)
summary(lm) # no subcascale in itself associated with RRST

Call:
lm(formula = metascore ~ group + ctq_physab + ctq_emoab + ctq_sexab + 
    ctq_emoneg + ctq_physneg, data = df)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.231470 -0.027816 -0.002843  0.039078  0.143323 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.5691398  0.0163848  34.736   <2e-16 ***
group1      -0.0132895  0.0122716  -1.083    0.282    
ctq_physab  -0.0002136  0.0024535  -0.087    0.931    
ctq_emoab    0.0028306  0.0021170   1.337    0.185    
ctq_sexab    0.0019599  0.0017353   1.129    0.262    
ctq_emoneg  -0.0024935  0.0017865  -1.396    0.166    
ctq_physneg -0.0004406  0.0023150  -0.190    0.850    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.05449 on 84 degrees of freedom
Multiple R-squared:  0.07209,   Adjusted R-squared:  0.005815 
F-statistic: 1.088 on 6 and 84 DF,  p-value: 0.3765
# predciting  traitPE uing RRST and MAIA with covariates
lm<-lm(traitPE_maiaRRST ~ group + sex + psychotropic_medication + anx_dep_SUM + ctq_total + sdq,
  data = df)
summary(lm) #  SEX, MEDICATION and ANX-DEP sum score sig associated with traitPE-maia

Call:
lm(formula = traitPE_maiaRRST ~ group + sex + psychotropic_medication + 
    anx_dep_SUM + ctq_total + sdq, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.71628 -0.73925  0.07294  0.80040  2.78511 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)               1.118175   0.466282   2.398  0.01875 *  
group1                   -0.204938   0.311930  -0.657  0.51302    
sex2                     -0.521139   0.271396  -1.920  0.05831 .  
psychotropic_medication1  1.002053   0.357972   2.799  0.00638 ** 
anx_dep_SUM              -0.030294   0.007368  -4.112 9.28e-05 ***
ctq_total                 0.002098   0.007180   0.292  0.77091    
sdq                       0.012108   0.012476   0.970  0.33467    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.107 on 82 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.2209,    Adjusted R-squared:  0.1639 
F-statistic: 3.876 on 6 and 82 DF,  p-value: 0.001867
# predciting  traitPE uing RRST and MAIA with covariates
lm<-lm(traitPE_iasRRST ~ group + sex + psychotropic_medication + anx_dep_SUM + ctq_total + sdq,
  data = df)
summary(lm) #  MEDICATION sig associated with traitPE-ias

Call:
lm(formula = traitPE_iasRRST ~ group + sex + psychotropic_medication + 
    anx_dep_SUM + ctq_total + sdq, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.1760 -0.7174  0.1740  0.6745  2.7776 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)  
(Intercept)               0.6802261  0.4918145   1.383   0.1704  
group1                   -0.1353882  0.3290103  -0.412   0.6818  
sex2                     -0.4335953  0.2862573  -1.515   0.1337  
psychotropic_medication1  0.7977286  0.3775735   2.113   0.0377 *
anx_dep_SUM              -0.0058473  0.0077715  -0.752   0.4540  
ctq_total                -0.0005382  0.0075729  -0.071   0.9435  
sdq                      -0.0112170  0.0131594  -0.852   0.3965  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.168 on 82 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.08506,   Adjusted R-squared:  0.01811 
F-statistic:  1.27 on 6 and 82 DF,  p-value: 0.28

4) Symptom Severity (FND only)

# Compute the correlation matrix with variables of symptom severity to see whether they correlate
# Select the desired variables
severity_variables <- df_FND[, c("sfmdrs", "cgi", "sdq", "sum_positive_signs", "sum_severity")]

# Compute the correlation matrix with rcorr
correlation_results <- Hmisc::rcorr(as.matrix(severity_variables))

# Extract correlation coefficients
cor_matrix <- correlation_results$r

# Extract p-values
p_matrix <- correlation_results$P

# Print results
cor_matrix
                       sfmdrs       cgi        sdq sum_positive_signs
sfmdrs             1.00000000 0.5841928 0.05625074          0.1716090
cgi                0.58419281 1.0000000 0.27890677          0.5531214
sdq                0.05625074 0.2789068 1.00000000          0.4738588
sum_positive_signs 0.17160897 0.5531214 0.47385879          1.0000000
sum_severity       0.17270871 0.4374258 0.06696646          0.1482822
                   sum_severity
sfmdrs               0.17270871
cgi                  0.43742577
sdq                  0.06696646
sum_positive_signs   0.14828221
sum_severity         1.00000000
p_matrix
                         sfmdrs          cgi         sdq sum_positive_signs
sfmdrs                       NA 3.136863e-05 0.716850710       2.653468e-01
cgi                3.136863e-05           NA 0.066747898       9.864238e-05
sdq                7.168507e-01 6.674790e-02          NA       1.157504e-03
sum_positive_signs 2.653468e-01 9.864238e-05 0.001157504                 NA
sum_severity       2.622499e-01 2.985648e-03 0.665812099       3.367527e-01
                   sum_severity
sfmdrs              0.262249862
cgi                 0.002985648
sdq                 0.665812099
sum_positive_signs  0.336752703
sum_severity                 NA
# SIG correlations of different symptom severity scores
### CGI & SFMDRS corr: r= 0.5841928 , p=3.136863e-05
### CGI & SUM_POS corr: r= 0.5531214  2, p=9.864238e-05
### SDQ & SUM_POS corr: r=0.4738588, p=1.157504e-03
### CGI & SUM_SEVERITY corr: r=0.43742577 , p=0.002985648

#SFMDRS
sfmdrslm <- lm(formula = sfmdrs ~  rrst_sensitivity + traitPE_iasRRST + traitPE_maiaRRST + metascore,
               data = df_FND) 
summary(sfmdrslm)   # no sig association

Call:
lm(formula = sfmdrs ~ rrst_sensitivity + traitPE_iasRRST + traitPE_maiaRRST + 
    metascore, data = df_FND)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.5557 -4.4430 -0.7832  2.8633 17.5487 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)  
(Intercept)       15.15210    8.12003   1.866   0.0696 .
rrst_sensitivity   0.98623    0.86030   1.146   0.2586  
traitPE_iasRRST    0.07545    0.83708   0.090   0.9286  
traitPE_maiaRRST   1.22722    0.96790   1.268   0.2123  
metascore        -20.23952   15.74488  -1.285   0.2062  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.149 on 39 degrees of freedom
Multiple R-squared:  0.06056,   Adjusted R-squared:  -0.03579 
F-statistic: 0.6285 on 4 and 39 DF,  p-value: 0.645
#CGI
cgilm <- lm(formula = cgi ~ rrst_sensitivity + traitPE_iasRRST + traitPE_maiaRRST + metascore,
            data = df_FND) 
summary(cgilm) # no sig assocaition

Call:
lm(formula = cgi ~ rrst_sensitivity + traitPE_iasRRST + traitPE_maiaRRST + 
    metascore, data = df_FND)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.4075 -1.2088 -0.1394  0.9222  2.9767 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)  
(Intercept)        5.0955     2.0335   2.506   0.0165 *
rrst_sensitivity   0.1073     0.2154   0.498   0.6213  
traitPE_iasRRST   -0.3201     0.2096  -1.527   0.1348  
traitPE_maiaRRST   0.4347     0.2424   1.793   0.0807 .
metascore         -4.8523     3.9431  -1.231   0.2258  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.54 on 39 degrees of freedom
Multiple R-squared:  0.1309,    Adjusted R-squared:  0.0418 
F-statistic: 1.469 on 4 and 39 DF,  p-value: 0.2302
#Sum score pos signs  (pos signs like wekaness in deltoid/stenocleido etc, co-contractions, give way weakness, hoover, hyposensitivity)
sum_pos_lm <- lm(formula = sum_positive_signs ~  rrst_sensitivity + metascore + ias + maia_total,
                 data = df_FND) 
summary(sum_pos_lm)   # p=0.0556 for rrst_sensitivity

Call:
lm(formula = sum_positive_signs ~ rrst_sensitivity + metascore + 
    ias + maia_total, data = df_FND)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.7530  -3.6712   0.5005   3.0771  25.9158 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)  
(Intercept)      11.79851    9.80874   1.203   0.2363  
rrst_sensitivity -1.30567    0.66170  -1.973   0.0556 .
metascore        10.72469   19.03536   0.563   0.5764  
ias              -0.04154    0.06807  -0.610   0.5452  
maia_total        0.31711    0.20353   1.558   0.1273  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.435 on 39 degrees of freedom
Multiple R-squared:  0.1428,    Adjusted R-squared:  0.05489 
F-statistic: 1.624 on 4 and 39 DF,  p-value: 0.1875
#Sum score severity (motor weakness and seizure frequency)
sum_pos_lm <- lm(formula = sum_severity ~  rrst_sensitivity + metascore + ias + maia_total,
                 data = df_FND) 
summary(sum_pos_lm)  # p= 0.069 for metascore

Call:
lm(formula = sum_severity ~ rrst_sensitivity + metascore + ias + 
    maia_total, data = df_FND)

Residuals:
    Min      1Q  Median      3Q     Max 
-52.138 -15.445  -6.742   5.041 245.056 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)  
(Intercept)       100.09921   56.43243   1.774   0.0839 .
rrst_sensitivity   -1.03365    3.80692  -0.272   0.7874  
metascore        -204.72797  109.51575  -1.869   0.0691 .
ias                 0.08866    0.39160   0.226   0.8221  
maia_total          1.42149    1.17099   1.214   0.2321  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 42.77 on 39 degrees of freedom
Multiple R-squared:  0.09938,   Adjusted R-squared:  0.007003 
F-statistic: 1.076 on 4 and 39 DF,  p-value: 0.3816

5) Visualization of sig association

df_plot <- df %>%
  mutate(sex = factor(sex, levels = c("1", "2"), labels = c("female", "male"))) %>%
  mutate(group = factor(group, levels = c("0", "1"), labels = c("HC", "FND")))

# SDQ-RRST_sensitivity
RRST_sdq <- ggplot(df_plot, aes(x = sdq, y = rrst_sensitivity, color = as.factor(group))) +
  geom_point() +  # Scatterplot for both groups
  geom_smooth(method = "lm", se = TRUE) +  # Add regression lines and SE
  labs(x = "Dissociative Score (SDQ-20)", 
       y = "Respiratory Sensitivity (RRST)", 
       title = "Respiratory Interoception (FND vs HC)",
       color = "Group") +  # Label for the legend
  theme_minimal() +
  scale_color_manual(values = c("HC" = "#868686FF", "FND" = "#E95248")) + #custom colors
  theme(panel.grid = element_blank())  # Remove grid lines

RRST_sdq
`geom_smooth()` using formula = 'y ~ x'

# CTQ-RRST_sensitivity
RRST_ctq <- ggplot(df_plot, aes(x = ctq_total, y = rrst_sensitivity, color = as.factor(group))) +
  geom_point() +  # Scatterplot for both groups
  geom_smooth(method = "lm", se = TRUE) +  # Add regression lines and SE
  labs(x = "Experience of Childhood Trauma (CTQ)", 
       y = "Respiratory Sensitivity (RRST)", 
       title = "Respiratory Interoception (FND vs HC)",
       color = "Group") +  # Label for the legend
  theme_minimal() +
  scale_color_manual(values = c("HC" = "#868686FF", "FND" = "#E95248")) + #custom colors
  theme(panel.grid = element_blank())  # Remove grid lines
RRST_ctq
`geom_smooth()` using formula = 'y ~ x'

# VISUALIZATION PER CTQ SUBSCORE that was foudn sig between groups
# Physical_Abuse-RRST_sensitivity
RRST_physab <- ggplot(df_plot, aes(x = ctq_physab, y = rrst_sensitivity, color = as.factor(group))) +
  geom_point() +  # Scatterplot for both groups
  geom_smooth(method = "lm", se = TRUE) +  # Add regression lines and SE
  labs(x = "Physical Abuse in Childhood (CTQ subscale)", 
       y = "Respiratory Sensitivity (RRST)", 
       title = "Respiratory Interoception (FND vs HC)",
       color = "Group") +  # Label for the legend
  theme_minimal() +
  scale_color_manual(values = c("HC" = "#868686FF", "FND" = "#E95248")) + #custom colors
  theme(panel.grid = element_blank())  # Remove grid lines
RRST_physab
`geom_smooth()` using formula = 'y ~ x'

# Emotional_Neglect-RRST_sensitivity
RRST_emoneg <- ggplot(df_plot, aes(x = ctq_emoneg, y = rrst_sensitivity, color = as.factor(group))) +
  geom_point() +  # Scatterplot for both groups
  geom_smooth(method = "lm", se = TRUE) +  # Add regression lines and SE
  labs(x = "Emotional Neglegt in Childhood (CTQ subscale)", 
       y = "Respiratory Sensitivity (RRST)", 
       title = "Respiratory Interoception (FND vs HC)",
       color = "Group") +  # Label for the legend
  theme_minimal() +
  scale_color_manual(values = c("HC" = "#868686FF", "FND" = "#E95248")) + #custom colors
  theme(panel.grid = element_blank())  # Remove grid lines
RRST_emoneg
`geom_smooth()` using formula = 'y ~ x'

# Emotional_Abuse-RRST_sensitivity
RRST_emoab <- ggplot(df_plot, aes(x = ctq_emoab, y = rrst_sensitivity, color = as.factor(group))) +
  geom_point() +  # Scatterplot for both groups
  geom_smooth(method = "lm", se = TRUE) +  # Add regression lines and SE
  labs(x = "Emotional Abuse in Childhood (CTQ subscale)", 
       y = "Respiratory Sensitivity (RRST)", 
       title = "Respiratory Interoception (FND vs HC)",
       color = "Group") +  # Label for the legend
  theme_minimal() +
  scale_color_manual(values = c("HC" = "#868686FF", "FND" = "#E95248")) + #custom colors
  theme(panel.grid = element_blank())  # Remove grid lines
RRST_emoab
`geom_smooth()` using formula = 'y ~ x'

# Sex-RRST_sensitivity
RRST_sex <- ggplot(df_plot, aes(x = as.factor(sex), y = rrst_sensitivity, fill = as.factor(group))) +
  geom_boxplot() +  # Create the box plot
  geom_jitter(position = position_jitterdodge(jitter.width = 0.5, dodge.width = 0.75), color = "black", size = 1.5, alpha = 0.6) +  # Add jittered points in black
  labs(x = "Sex", 
       y = "TraitPE (HBT-IAS)", 
       title = "Trait-Prediction-Error HBT-IAS (FND vs HC)",
       fill = "Group") +  # Label for the legend
  theme_minimal() +
  scale_fill_manual(values = c("HC" = "#868686FF", "FND" = "#E95248"))+ #custom colors
  theme(panel.grid = element_blank())  # Remove grid lines
RRST_sex

# SDQ-metascore
meta_sdq <- ggplot(df_plot, aes(x = sdq, y = metascore, color = as.factor(group))) +
  geom_point() +  # Scatterplot for both groups
  geom_smooth(method = "lm", se = TRUE) +  # Add regression lines and SE
  labs(x = "Dissociative Score (SDQ-20)", 
       y = "Respiratory Metacognition (AUC)", 
       title = "Respiratory Interoception (FND vs HC)",
       color = "Group") +  # Label for the legend
  theme_minimal() +
  scale_color_manual(values = c("HC" = "#868686FF", "FND" = "#E95248")) + #custom colors
  theme(panel.grid = element_blank())  # Remove grid lines
meta_sdq
`geom_smooth()` using formula = 'y ~ x'

# CTQ-metascore
meta_ctq <- ggplot(df_plot, aes(x = ctq_total, y = metascore, color = as.factor(group))) +
  geom_point() +  # Scatterplot for both groups
  geom_smooth(method = "lm", se = TRUE) +  # Add regression lines and SE
  labs(x = "Experience of Childhood Trauma (CTQ)", 
       y = "Respiratory Metacognition (AUC)", 
       title = "Respiratory Interoception (FND vs HC)",
       color = "Group") +  # Label for the legend
  theme_minimal() +
  scale_color_manual(values = c("HC" = "#868686FF", "FND" = "#E95248"))+ #custom colors
  theme(panel.grid = element_blank())  # Remove grid lines
meta_ctq
`geom_smooth()` using formula = 'y ~ x'

# traitPEmaia-Anx-Dep
traitPEmaia_anxdep <- ggplot(df_plot, aes(x = anx_dep_SUM, y = traitPE_maiaRRST, color = as.factor(group))) +
  geom_point() +  # Scatterplot for both groups
  geom_smooth(method = "lm", se = TRUE) +  # Add regression lines and SE
  labs(x = "Anxiety & Depression Sum Score (STAI-T & BDI-II)", 
       y = "Trait-PredictionError (RRST-MAIA)", 
       title = "Interoceptive Trait-Prediction Error (FND vs HC)",
       color = "Group") +  # Label for the legend
  theme_minimal() +
  scale_color_manual(values = c("HC" = "#868686FF", "FND" = "#E95248")) + #custom colors
  theme(panel.grid = element_blank())  # Remove grid lines
traitPEmaia_anxdep
`geom_smooth()` using formula = 'y ~ x'

Take Home Message

After excluding one outlier (FND patient > 2.5 SD from mean, and outlier in terms of value, but also in testing the model fit in terms of normality/homoscedasticity of resiudaly) we have N=91.

Baseline: FND patients have higher intake of psychotropic medication (p<0.001), higher BMI (p=0.019), higher Depression (p=4.2e-09) and Anxiety Scores (State: p=4.9e-06, Trait: p=4.1e-04), higher Experience of Childhood Trauma (p=0.0014) and Dissociation (p=1.4e-10) compared to controls.

Interoception: FND patients show to have lower interoceptive sensibility, juding from MAIA questionnaire (p=4.4e-04) and IAS (p=0.018) compared to controls The traitPE of RRST and these two questionnaires did not differ sig between the groups.

FND patients have lower interoceptive sensitivity (p=0.032), but no difference in metacognition (p=0.9), compared to controls.

SDQ is associated with GROUP, when controlling for clinical and interoceptive variables (p=0.00147).

Clinical Covariates to include in the model were: Group, Sex, Intake of psychotropic_medication (yes/no), Anxiety/Depression Sum Score, CTQ total score, SDQ-20 score: Male sex (t=2.456, p=0.01616), no intake of psychotropic medication (-2.667, p=0.00921), experience of more childhood trauma (t=2.277, p=0.02541) and lower dissociation (t=-2.894, p=0.00487) is associated with higher RRST sensitivity, when controlling for all other variables. Lower dissociation (t=-2.747, p=0.00739) and with weak evidence also Childhood Trauma (t=1.934, p=0.05651) were associated with higher metascores, when controlling for all other variables. NOTE: the CTQ subscores would not show visually / nor statistically a specific contirbution to the RRST sensitivity or metacognition. The intake of psychotropic medication (t=2.799, p=0.00638) and lower scores in Anxiety/Depression (t=-4.112, p=9.28e-05), and with weak evidence also female sex (t=-1.920, p=0.05831) are associated with higher TraitPrediction Error using RRST-sensitivity and the MAIA questionnaire, when controlling all other variables. The intake of psychotropic medication (t=2.113, p=0.0377) is associated with higher TraitPrediction Error using RRST-sensitivity and the IAS questionnaire, when controlling all other variables.

Symptom Severity There is weak evidence, that lower RRST-sensitivity is associated with more positive signs in FND (t=-1.973, p=0.0556).

From the visualization it is apparent that the RRST-Sensitivity & Metacognition in Controls is not dependent on dissociation scores, while for patients they are lower for higher dissociation scores (but note also the larger spread for patients only). While there is no linear association for CTQ and Metacognition, it iseems that for FND only, the more experiences of Childhood Trauma, the higher the metacognitive score. Finally, the Trait-Prediction Error (RRST-MAIA) shows that controls seem to overestimate themselve sin terms of interoceptive performance when scoring low on Anxiety&Depression, while they start to underestimate their performance with higher Anxiety&Depression scores. For FND, this trend is definitely smaller; independent of affective smyptoms, FND patients seem to be closer to a zero in Trait-Prediction Error (using RRST performance and MAIA self report).